#include "stdafx.h"
#include "RayTracer.h"
#include <cfloat>
#include <cmath>
#include <cstdlib>
#include "tinyxml.h"
using namespace std;

namespace raytracer
{
	typedef boost::ptr_vector<Light> LightV;
	typedef boost::ptr_vector<Primitive> PrimitiveV;
	typedef unsigned char uchar;

	//a very little number 
	const float Epsilon = 0.001f;
	//This version just ignores the filename and generate scene with code
	bool Scene::InitScene(const string& filename)
	{
		//clean
		delete [] mColorBuffer;
		mLights.clear();
		mPrimitives.clear();

		// load scene from file

		TiXmlDocument doc(filename.c_str());
		if(doc.LoadFile())
		{
			const TiXmlElement* rootEle = doc.RootElement();
			mWX = (float)atof(rootEle->Attribute("WorldX"));
			mWY = (float)atof(rootEle->Attribute("WorldY"));
			mWWidth = (float)atof(rootEle->Attribute("WorldWidth"));
			mWHeight = (float)atof(rootEle->Attribute("WorldHeight"));
			mWidth = atoi(rootEle->Attribute("ImageWidth"));
			mHeight = atoi(rootEle->Attribute("ImageHeight"));

			//set eye position
			const TiXmlElement* eyePosEle = rootEle->FirstChildElement("EyePosition");
			assert(eyePosEle);
			mEyePos = vector3(float(atof(eyePosEle->Attribute("x"))),
				float(atof(eyePosEle->Attribute("y"))),
				float(atof(eyePosEle->Attribute("z"))));

			//set Ambient
			const TiXmlElement* ambEle = rootEle->FirstChildElement("Ambient");
			assert(ambEle);
			mAmbient = vector3(float(atof(ambEle->Attribute("r"))),
				float(atof(ambEle->Attribute("g"))),
				float(atof(ambEle->Attribute("b"))));

			//set attenuation
			const TiXmlElement* attEle = rootEle->FirstChildElement("Attenuation");
			assert(attEle);
			mAttIdx = vector3(float(atof(attEle->Attribute("a"))),
				float(atof(attEle->Attribute("b"))),
				float(atof(attEle->Attribute("c"))));

			//other configs
			mMaxStep = 6;
			mColorBuffer = new int[mWidth*mHeight];

			const TiXmlElement* ele = rootEle->FirstChildElement();
			while(ele)
			{

				//light
				if(!strcmp("Light",ele->Value()))
				{
					const TiXmlElement* lightEle = ele;
					const TiXmlElement* colorEle,*centerEle;
					vector3 color,center;
					Light* pLight;
					colorEle = lightEle->FirstChildElement("Color");
					centerEle = lightEle->FirstChildElement("Center");
					assert(colorEle && centerEle);
					color = vector3(float(atof(colorEle->Attribute("r"))),
						float(atof(colorEle->Attribute("g"))),
						float(atof(colorEle->Attribute("b"))));
					center = vector3(float(atof(centerEle->Attribute("x"))),
						float(atof(centerEle->Attribute("y"))),
						float(atof(centerEle->Attribute("z"))));
					pLight = new Light(color,Light::Point);
					pLight->SetCenter(center);
					mLights.push_back(pLight);
				}
				else if(!strcmp("Sphere",ele->Value()))
				{
					mPrimitives.push_back(new Sphere(ele));
				}
				else if(!strcmp("Plane",ele->Value()))
				{
					mPrimitives.push_back(new Plane(ele));
				}
				else if(!strcmp("Polyhedron",ele->Value()))
				{
					Material m;
					m.Set(ele->FirstChildElement("Material"));
					const TiXmlElement* triListEle = ele->FirstChildElement("Triangles");
					assert(triListEle);
					const TiXmlElement* triEle = triListEle->FirstChildElement("Triangle");
					while(triEle)
					{
						Triangle* pTri = new Triangle(triEle);
						pTri->SetMaterial(m);
						mPrimitives.push_back(pTri);
						triEle = triEle->NextSiblingElement("Triangle");
					}
				}
				else if(!strcmp("Cylindar",ele->Value()))
				{
					mPrimitives.push_back(new Cylindar(ele));
				}
				ele = ele->NextSiblingElement();
			}
			return true;
		}
		else
			return false;



		//mWidth=800;
		//mHeight=600;
		//mWX=-4,mWY=3;
		//mWWidth=8,mWHeight=6;
		//mMaxStep = 6;
		//mColorBuffer = new int[mHeight*mWidth];
		//mEyePos=vector3(0,0,-5);
		//mAmbient=vector3(0,0,0);
		//mAttIdx=vector3(0,0,1);




		////add lights
		//Light* light = new Light(vector3( 0.4f, 0.4f, 0.4f),Light::Point);
		//light->SetCenter(vector3( 0, 5, 5 ));
		//mLights.push_back(light);
		//light = new Light(vector3(0.6f, 0.6f, 0.8f),Light::Point);
		//light->SetCenter(vector3( -3, 5, 1 ));
		//mLights.push_back(light);

		////add primitives
		//Material m;
		//m.SetAmbientReflection(0.5f);
		//m.SetDiffuseReflection(1.f);
		//m.SetSpecularPower(20);
		//m.SetSpecular(1.f);
		////big sphere
		//m.SetDiffuseReflection(0.8f);
		//m.SetDiffuse(vector3(0.7f, 0.7f, 1.f ));
		//m.SetReflection(0.2f);
		//m.SetRefractionIndex(1.3f);
		//m.SetRefraction(0.8f);
		//m.SetSpecular(1.f);
		//mPrimitives.push_back(new Sphere(m,vector3( 2, 0.8f, 3 ), 2.5f ));
		////small sphere
		//m.SetDiffuse(vector3( 0.7f, 0.7f, 1.0f ));
		//m.SetReflection(0.5f);
		//m.SetRefraction(0.f);
		//m.SetRefractionIndex(1.3f);
		//m.SetDiffuseReflection(0.3f);
		//mPrimitives.push_back(new Sphere(m,vector3( -5.5f, -0.5f, 7  ), 2 ));
		////behind sphere
		//m.SetDiffuse(vector3(0.8f,0.2f,0.2f));
		//m.SetDiffuseReflection(1.f);
		//m.SetReflection(0.f);
		//m.SetRefraction(0.f);
		//m.SetSpecular(1.f);
		//mPrimitives.push_back(new Sphere(m,vector3(2,0.8f,7),0.3f));
		////ground
		//m.SetDiffuse(vector3(0.4f, 0.3f, 0.3f));
		//m.SetReflection(0);
		//m.SetRefraction(0);
		//m.SetSpecular(0);
		//m.SetDiffuseReflection(1.f);
		//mPrimitives.push_back(new Plane(m,vector3( 0, 1, 0 ), 4.4f));
		////a triangle
		//vector3 triPoints[3] = {vector3(-3,0,0),vector3(-2,0,0),vector3(-2.5f,-1,0)};
		//m.SetDiffuse(vector3(0.2f,0.2f,0.9f));
		//m.SetDiffuseReflection(1.f);
		//m.SetReflection(0.f);
		//m.SetRefraction(0.f);
		//mPrimitives.push_back(new Triangle(m,triPoints));
	}

	//Render the scene
	void Scene::Render()
	{

		/*long curIndex = 0;
		for(int j=0;j<mHeight;j++)
		for(int i=0;i<mWidth;i++,curIndex++)
		mColorBuffer[curIndex]=0x00FFFF;*/

		long curBufferIdx = 0;
		float deltaX = mWWidth/mWidth,deltaY = mWHeight/mHeight;
		//current x,y in world
		double wx,wy=mWY-mWHeight;
		for(int h=0;h<mHeight;h++,wy+=deltaY)
		{
			wx = mWX;
			for(int w=0;w<mWidth;w++,curBufferIdx++,wx += deltaX)
			{
				//ray from position
				vector3 scrPos(float(wx),float(wy),0);
				Ray srcRay(mEyePos/*scrPos*/,scrPos-mEyePos);
				vector3 color;	
				float dist;
				if((w==500) && (h==240))
				{
					int a = 0;
					a++;
				}
				RayTrace(srcRay,color,1,1,dist);

				//this version has no reflection just plain view
				//int prIdx = -1;
				//double minDistance = DBL_MAX;	
				////find the nearest primitive to intersect with the ray
				//for(PrimitiveV::size_type i=0;i<mPrimitives.size();i++)
				//{
				//double dis = mPrimitives[i].Intersect(srcRay);
				//if((dis > Epsilon) && (dis<minDistance))
				//{
				//prIdx=int(i);
				//minDistance = dis;
				//}
				//}

				//if(prIdx != -1)
				//{
				//const Material& m = mPrimitives[prIdx].GetMaterial();
				//color.r = m.GetDiffuse().r;
				//color.g = m.GetDiffuse().g;
				//color.b = m.GetDiffuse().b;
				//}				

				mColorBuffer[curBufferIdx] = uchar(min(color.r*255,255)) + ((uchar(min(color.g*255,255)))<<8) + ((uchar(min(color.b*255,255)))<<16);
			}
		}
	}

	void Sphere::Intersect(const Ray& ray,bool* pIsIntted,double* pLen,vector3* pInttPt,vector3* pNormal)const
	{
		assert(pIsIntted);
		vector3 dist = mCenter-ray.GetOrigin();
		float B = dist.Dot(ray.GetDirection());
		float D = B*B - dist.SqrLength() + mRadius*mRadius;


		if(D<0)
		{
			*pIsIntted=false;
			return;
		}

		double len;
		vector3 inttPt,normal;

		float t0=B - sqrtf(D),t1=B+sqrtf(D);
		if(t0>Epsilon)
		{
			*pIsIntted=true;
			len=t0;			
		}
		else if(t1 >Epsilon)
		{
			*pIsIntted=true;
			len=t1;
		}
		else
		{
			*pIsIntted=false;
			return;
		}


		if(pLen)
			*pLen=len;

		if(!pInttPt && !pNormal)
			return;

		inttPt=(ray.GetOrigin()+float(len)*ray.GetDirection());
		//since it may intersect from the inside
		normal=(inttPt-mCenter);
		if(normal.Dot(ray.GetDirection()) >= Epsilon)
			normal *= -1.f;
		if(pInttPt)
			*pInttPt=inttPt;		
		if(pNormal)
		{
			*pNormal=normal;
			pNormal->Normalize();
		}

	}
	/*double Sphere::Intersect(const Ray& ray)const
	{
	vector3 dist = mCenter-ray.GetOrigin();
	vector3 rayDir = ray.GetDirection();
	rayDir.Normalize();

	float B = dist.Dot(rayDir);
	if(B <= 0)
	return -1;

	float D = B*B - dist.SqrLength() + mRadius*mRadius;
	if(D < 0)
	return -1;

	float t0 = B - sqrtf(D);
	if(t0>Epsilon)
	return t0;

	float t1 = B + sqrtf(D);
	if(t1>Epsilon)
	return t1;

	return -1;
	}*/

	/*
	An recursive implementation of ray trace
	depth is the depth of the recursion
	refrIndex is the refraction index
	dist is the distance the ray travels, which will be taken into account when using transparency
	*/
	Primitive* Scene::RayTrace(Ray& ray,vector3& color, int depth, float refrIndex, float& dist)
	{
		//clear values
		color = vector3();
		dist=0;

		//the index of the nearest primitive
		int prIdx = -1;
		double minDistance = DBL_MAX;	
		//intersection point and normal
		vector3 nPoint,priNormal;
		//ray.NormalizeDirection();
		//find the nearest primitive to intersect with the ray
		for(PrimitiveV::size_type i=0;i<mPrimitives.size();i++)
		{
			bool isIntted;
			double dis;
			vector3 inttPt,inttNrm;
			mPrimitives[i].Intersect(ray,&isIntted,&dis,&inttPt,&inttNrm);
			if(isIntted && (dis<minDistance))
			{
				prIdx=int(i);
				minDistance = dis;
				nPoint=inttPt;
				priNormal=inttNrm;
			}
		}

		if(prIdx == -1)
		{
			//if it can not hit a object,we test wheter it can hit a ball centered at light
			/*Sphere s(Material(),vector3(),1.f);;
			bool isIntted;
			for(LightV::size_type i=0;i<mLights.size();i++)
			{
			s.SetCenter(mLights[i].GetCenter());
			s.Intersect(ray,&isIntted);
			if(isIntted)
			{
			color=vector3(1,1,1);
			break;
			}
			}*/
			return NULL;
		}
		const Material& material = mPrimitives[prIdx].GetMaterial();
		//ambient light
		color = (mAmbient * material.GetAmbientReflection()) * material.GetDiffuse();

		//the intersection point
		//vector3 nPoint = ray.GetOrigin() + float(minDistance)*ray.GetDirection();
		//primitive normal at intersection point
		//vector3 priNormal = mPrimitives[prIdx].GetNormalAt(nPoint);


		//calculate contributions of lights to the intersection point
		for(LightV::size_type i = 0;i<mLights.size();i++)
		{
			//light ray at intersection point
			Ray lRay;
			//the length of lien segment between light source and point
			double rayLen = DBL_MAX;
			//to add other light type support
			if(mLights[i].GetLightType()==Light::Point)
			{
				lRay.SetOrigin(mLights[i].GetCenter());
				lRay.SetDirection(mLights[i].GetDirAt(nPoint));
				rayLen=(mLights[i].GetCenter()-nPoint).Length();
				/*if((priNormal.Dot(lRay.GetDirection())>Epsilon) || (rayLen <= Epsilon))
				continue;*/
			}




			Ray shadowTestRay(nPoint,lRay.GetDirection()*-1);
			//test whether the light is shaded by an primitive
			bool isInShadow = false;
			for(PrimitiveV::size_type j = 0;(j < mPrimitives.size());j++)
			{
				bool isIntted;
				double len;
				mPrimitives[j].Intersect(shadowTestRay,&isIntted,&len);
				if(isIntted && (len < rayLen))
				{
					isInShadow = true;
					break;
				}
			}
			if(isInShadow)
				continue;


			//attenuation
			float atten = 1.f;
			//for point light the ray is a line segment			
			if(Light::Point==mLights[i].GetLightType())
			{
				atten *= (float)min(1,1/(mAttIdx.cell[0]*rayLen*rayLen+mAttIdx.cell[1]*rayLen+mAttIdx.cell[2]));
				atten = abs(atten);
				//lambert *= abs(lRay.GetDirection().Dot(priNormal));
			}


			//diffuse part
			color += atten * mLights[i].GetColor()*material.GetDiffuseReflection()*material.GetDiffuse()*abs(lRay.GetDirection().Dot(priNormal));
			//specular part
			if(material.GetSpecular() > Epsilon)
			{
				vector3 lRefl = lRay.GetDirection() - 2 * lRay.GetDirection().Dot(priNormal) * priNormal;
				lRefl.Normalize();
				color += atten * mLights[i].GetColor()*material.GetSpecular()
					*pow(abs(lRefl.Dot(ray.GetDirection())),material.GetSpecularPower());
			}			
		}

		//if the recursive step exceeds the max step, stop
		if(depth >= mMaxStep)
			return NULL;

		//reflection
		if(material.GetReflection() > Epsilon)
		{
			vector3 reflcolor;
			float dist;
			Ray reflray(nPoint,ray.GetDirection()-2*ray.GetDirection().Dot(priNormal)*priNormal);
			//reflray.SetOrigin(nPoint + Epsilon * reflray.GetDirection());
			RayTrace(reflray,reflcolor,depth+1,refrIndex,dist);
			color += material.GetReflection() * reflcolor;
		}
		//refraction
		if(material.GetRefraction() > Epsilon)
		{
			//the way to calculate refraction ray is found here http://www.devmaster.net/articles/raytracing_series/part3.php
			float ratio = refrIndex/material.GetRefractionIndex();
			float cosI = priNormal.Dot(ray.GetDirection());
			float cosT2 = 1.f - ratio*ratio*(1-cosI*cosI);
			if(cosT2 > Epsilon)
			{
				vector3 refrDir=ratio * ray.GetDirection() - (ratio*cosI + sqrtf(cosT2))*priNormal;
				vector3 refrColor;
				float dist;
				RayTrace(Ray(nPoint,refrDir),refrColor,depth+1,material.GetRefractionIndex(),dist);
				color += refrColor;
			}
		}

		/*color.r = min(1.f,color.r);
		color.g = min(1.f,color.g);
		color.b = min(1.f,color.b);*/
		return &(mPrimitives[prIdx]);
	}

	void Plane::Intersect(const Ray& ray,bool* pIsIntted,double* pLen,vector3* pInttPt,vector3* pNormal)const
	{
		assert(pIsIntted);

		double dot = mNormal.Dot(ray.GetDirection());
		if(abs(dot) > Epsilon)
		{
			double dist = -(mNormal.Dot(ray.GetOrigin()) + mDist)/dot;
			if(dist > Epsilon)//intersected
			{
				*pIsIntted=true;
				if(pLen)
					*pLen=dist;
				if(pInttPt)
					*pInttPt=ray.GetOrigin()+ float(dist)* ray.GetDirection();
				if(pNormal)
					*pNormal=(ray.GetDirection().Dot(mNormal)>=0)?(-1*mNormal):mNormal;
			}
			else
				*pIsIntted=false;

		}
		else//no intersection
			*pIsIntted=false;

	}

	//double Plane::Intersect(const Ray& ray)const
	//{
	//	double dot = mNormal.Dot(ray.GetDirection());
	//	if(abs(dot) > Epsilon)
	//	{
	//		double dist = -(mNormal.Dot(ray.GetOrigin()) + mDist)/dot;
	//		if(dist > Epsilon)
	//			return dist;
	//		else
	//			return -1;

	//	}
	//	else//no intersection
	//		return -1;
	//}

	//Polyhedron::Polyhedron(const Material&m,const std::vector<vector3>& points,const std::vector<std::vector<vector3>::size_type>& indexes)
	//	:Primitive(m),mPoints(points),mIndexes(indexes)
	//{
	//	//calculate normals
	//	for(IndexV::size_type i=0;i<mIndexes.size();i += 3)
	//	{
	//		mNormals.push_back((mPoints[mIndexes[i+1]]-mPoints[mIndexes[i]]).Cross(mPoints[mIndexes[i+2]]-mPoints[mIndexes[i]]));
	//	}		
	//}

	////The theory of the algorithm please refer to http://hi.baidu.com/chuckhao/blog/item/863e08a9cb658bf91e17a2a0.html
	//void Polyhedron::Intersect(const Ray& ray,bool* pIsIntted,double* pLen=NULL,vector3* pInttPt=NULL,vector3* pNormal=NULL)const
	//{
	//	assert(pIsIntted);
	//	for(IndexV::size_type i=0;i<mIndexes.size();i+=3)
	//	{
	//		const vector3& A = mPoints[mIndexes[i]];
	//		const vector3& B = mPoints[mIndexes[i+1]];
	//		const vector3& C = mPoints[mIndexes[i+2]];			
	//		const vector3& R = ray.GetOrigin();
	//		const vector3& D = ray.GetDirection();
	//		double detA = -1 * (A-B).Cross(C-A).Dot(D);
	//		if(abs(detA)<Epsilon)
	//			continue;
	//		
	//		double q = -1 * (A-R).Cross(C-A).Dot(D)/detA;
	//	}
	//	
	//}

	Triangle::Triangle(const Material& m,const vector3 points[3]):Primitive(m)
	{
		/*for(int i=0;i<3;i++)
		mCoords[i]=points[i];*/
		memcpy(mCoords,points,sizeof(vector3)*3);
		InitNormal();
	}

	void Triangle::Intersect(const Ray& ray,bool* pIsIntted,double* pLen,vector3* pInttPt,vector3* pNormal)const
	{
		const vector3& A = mCoords[0];
		const vector3& B = mCoords[1];
		const vector3& C = mCoords[2];			
		const vector3& R = ray.GetOrigin();
		const vector3& D = ray.GetDirection();
		double detA = (C-A).Cross(A-B).Dot(D);
		if(abs(detA)<Epsilon)
		{
			*pIsIntted = false;
			return;
		}

		double q = (C-A).Cross(A-R).Dot(D)/detA;
		if( (q<=Epsilon) || (q > 1))
		{
			*pIsIntted = false;
			return;
		}


		double r = (R-A).Cross(A-B).Dot(D)/detA;
		if((r<=Epsilon) || (r>1))
		{
			*pIsIntted = false;
			return;
		}

		if((q+r)>1)
		{
			*pIsIntted = false;
			return;
		}

		double t = (C-A).Cross(A-B).Dot(A-R)/detA;
		if(t<Epsilon)
		{
			*pIsIntted = false;
			return;
		}

		*pIsIntted = true;

		if(pLen)
			*pLen=t;

		if(pInttPt)
			*pInttPt = ray.GetOrigin() + ray.GetDirection() * float(t);
		if(pNormal)
			*pNormal = (mNormal.Dot(ray.GetDirection())>Epsilon)?(-1.f * mNormal):mNormal;
	}

	void Material::Set(const TiXmlElement* ele)
	{
		assert(ele);
		mAmbnRef = float(atof(
			ele->FirstChildElement("AmbientReflection")->Attribute("value")));
		mDiffRef = float(atof(
			ele->FirstChildElement("DiffuseReflection")->Attribute("value")));
		mSpecPow = atoi(
			ele->FirstChildElement("SpecularPower")->Attribute("value"));
		mSpecRef = float(atof(
			ele->FirstChildElement("Specular")->Attribute("value")));
		const TiXmlElement* diffEle = ele->FirstChildElement("Diffuse");
		assert(diffEle);
		diffEle->QueryFloatAttribute("r",&mDiffuse.r);
		diffEle->QueryFloatAttribute("g",&mDiffuse.g);
		diffEle->QueryFloatAttribute("b",&mDiffuse.b);
		mReflIdx = float(atof(
			ele->FirstChildElement("Reflection")->Attribute("value")));
		mRefrRef = float(atof(
			ele->FirstChildElement("Refraction")->Attribute("value")));
		mRefrIdx = float(atof(
			ele->FirstChildElement("RefractionIndex")->Attribute("value")));
	}

	Plane::Plane(const TiXmlElement* ele)
	{
		assert(ele);
		mMaterial.Set(ele->FirstChildElement("Material"));
		const TiXmlElement* argsEle = ele->FirstChildElement("Args");
		assert(argsEle);
		argsEle->QueryFloatAttribute("x",&mNormal.x);
		argsEle->QueryFloatAttribute("y",&mNormal.y);
		argsEle->QueryFloatAttribute("z",&mNormal.z);
		argsEle->QueryFloatAttribute("d",&mDist);
		mNormal.Normalize();
	}

	Sphere::Sphere(const TiXmlElement* ele)
	{
		assert(ele);
		mMaterial.Set(ele->FirstChildElement("Material"));
		const TiXmlElement* centerEle = ele->FirstChildElement("Center");
		const TiXmlElement* radiusEle = ele->FirstChildElement("Radius");
		assert(centerEle&&radiusEle);
		centerEle->QueryFloatAttribute("x",&mCenter.x);
		centerEle->QueryFloatAttribute("y",&mCenter.y);
		centerEle->QueryFloatAttribute("z",&mCenter.z);
		radiusEle->QueryFloatAttribute("value",&mRadius);
	}

	//no material, only points
	Triangle::Triangle(const TiXmlElement* ele)
	{
		assert(ele);
		const TiXmlElement* ptEle = ele->FirstChildElement("Point");
		for(int i=0;i<3;i++,ptEle = ptEle->NextSiblingElement("Point"))
		{
			assert(ptEle);
			ptEle->QueryFloatAttribute("x",&(mCoords[i].x));
			ptEle->QueryFloatAttribute("y",&(mCoords[i].y));
			ptEle->QueryFloatAttribute("z",&(mCoords[i].z));
		}
		InitNormal();
	}

	Cylindar::Cylindar(const TiXmlElement* ele)
	{
		assert(ele);
		mMaterial.Set(ele->FirstChildElement("Material"));
		ele->FirstChildElement("Height")->QueryFloatAttribute("value",&mHeight);
		ele->FirstChildElement("Radius")->QueryFloatAttribute("value",&mRadius);
		const TiXmlElement* bcEle = ele->FirstChildElement("Bottom");
		assert(bcEle);
		bcEle->QueryFloatAttribute("x",&mBCenter.x);
		bcEle->QueryFloatAttribute("y",&mBCenter.y);
		bcEle->QueryFloatAttribute("z",&mBCenter.z);
		const TiXmlElement* dirEle = ele->FirstChildElement("Direction");
		assert(dirEle);
		dirEle->QueryFloatAttribute("x",&mVDir.x);
		dirEle->QueryFloatAttribute("y",&mVDir.y);
		dirEle->QueryFloatAttribute("z",&mVDir.z);
		mVDir.Normalize();
	}

	void Cylindar::Intersect(const Ray& ray,bool* pIsIntted,double* pLen,vector3* pInttPt,vector3* pNormal)const
	{
		const vector3& r = ray.GetOrigin();
		const vector3& d = ray.GetDirection();

		float len1=FLT_MAX,len2=FLT_MAX;
		bool isIntted1 = false,isIntted2 = false;
		//test intersection with bottom face
		float tmp1 = d.Dot(mVDir);
		if(abs(tmp1) > Epsilon)
		{
			len1 = (mBCenter-r).Dot(mVDir)/float(tmp1);
			isIntted1 = ((len1>Epsilon)&&((r+len1*d-mBCenter).SqrLength()<(mRadius*mRadius)));
				

			//top center
			vector3 tCen = mBCenter + mHeight*mVDir;
			len2 = (tCen-r).Dot(mVDir)/tmp1;
			isIntted2 = ((len2>Epsilon)&&((r+len2*d-tCen).SqrLength()<float(mRadius*mRadius)));

			if(isIntted1)
			{
				if(isIntted2)
				{
					len1 = min(len1,len2);
				}
			}
			else if(isIntted2)
			{
				isIntted1 = true;
				len1 = len2;
			}
		}



		//test intersection with body
		isIntted2 = false;
		float a = d.Cross(mVDir).SqrLength();
		float b = (r-mBCenter).Cross(mVDir).Dot(d.Cross(mVDir));
		float c = (r-mBCenter).Cross(mVDir).SqrLength()-mRadius*mRadius;
		float height;
		if(abs(a)<Epsilon)
		{
			len2 = -c/(2*b);
			height = (r+d*float(len2)).Dot(mVDir);
			if((len2>Epsilon)&&(height>Epsilon)&&(height<=mHeight))
				isIntted2 = true;
		}
		else
		{
			float det = b*b-a*c;
			if(det>Epsilon)
			{
				float t0 = (-b-sqrt(det))/a;
				height = (r+d*float(t0)-mBCenter).Dot(mVDir);
				if((t0 < Epsilon) || (height<Epsilon)||(height>mHeight))
				{
					t0=FLT_MAX;
				}

				float t1 = (-b+sqrt(det))/a;
				height = (r+d*float(t1)-mBCenter).Dot(mVDir);
				if((t1 < Epsilon) || (height<Epsilon)||(height>mHeight))
				{
					t1=FLT_MAX;
				}

				if(FLT_MAX==t0)
				{
					if(FLT_MAX>t1)
					{
						isIntted2=true;
						len2=t1;
					}
				}
				else
				{
					isIntted2=true;
					len2 = min(t1,t0);
				}
			}
		}

		*pIsIntted = isIntted1 || isIntted2;

		if(!(*pIsIntted))
			return;
		//indicates whether intersects bottom face
		bool isInttBt;

		double len;
		vector3 inttPt,inttNrm;
		if(isIntted1 && isIntted2)
		{
			isInttBt = len1<len2;
		}
		else
			isInttBt = isIntted1?true:false;

		if(!isInttBt)
		{
			len = len2;
			inttPt = r+float(len)*d;
			inttNrm = inttPt-mBCenter-(inttPt-mBCenter).Dot(mVDir)*mVDir;
			inttNrm *= inttNrm.Dot(d)>0?-1.f:1.f;
		}
		else
		{
			len = len1;
			inttPt = r+float(len)*d;
			inttNrm = (d.Dot(mVDir)<0)?mVDir:(-1*mVDir);
		}

		if(pLen)
			*pLen=len;
		if(pInttPt)
			*pInttPt=inttPt;
		if(pNormal)
		{
			*pNormal=inttNrm;
			pNormal->Normalize();
		}
	}
}
